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ABSTRACT 

The Modified Garabedian-McFadden (MGM) design procedure was incorporated into an 
existing 2D multigrid Navier-Stokes airfoil analysis method. The resulting design method is 
an iterative procedure based on a residual-correction algorithm and permits the automated design 
of airfoil sections with prescribed surface pressure distributions. The new design method, MG- 
MGM, is demonstrated for several different transonic pressure distributions obtained from both 
symmetric and cambered airfoil shapes. The airfoil profiles generated with the MG-MGM code 
are compared to the original configurations to assess the capabilities of the inverse design method. 

INTRODUCTION 

The aerodynamic design of aircraft components is often carried out by means of one of the 
following four approaches: a) cut-and-try analysis, b) indirect methods, c) optimization techniques, 
and d) inverse design techniques. Unlike the cut-and-try method, the latter three design techniques 
are far more automated, and can significantly reduce the overall engineering effort and calendar 
time required for developing aircraft components and configurations with improved aerodynamic 
performance or aerodynamic interference characteristics. 

A common design approach is to specify, a priori, surface pressure distributions that have fa- 
vorable aerodynamic characteristics at given freestream conditions. For example, an appropriately 
chosen pressure distribution can be used to achieve certain desired lift and moment coefficient 
goals, while a “weak-shock” or “shock-free” distribution can be used to minimize wave drag 
performance penalties. The automated design procedure is then used to generate, as efficiently 
as possible, the configuration geometry which will cause the specified pressures to exist on the 
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designed component. Obviously, the use of these automated design methods requires that the 
aerodynamicist can specify, a priori, the desired pressure distributions for a particular application. 

The most widely used aerodynamic design procedures for transonic-flow applications seem to 
be based upon potential-flow Computational Fluid Dynamics (CFD) methods. 1-5 This trend is most 
likely due to the relative low cost, in terms of computer-resource requirements, that is demonstrated 
by CFD methods based on the Transonic Small Disturbance (TSD) equation or the Full Potential 
equation (FPE). In the past decade, however, considerable interest has been demonstrated in the 
use of higher-order CFD methods such as the Euler equations and the Reynolds-averaged Navier- 
Stokes equations (RANS) for aerodynamic analyses in a variety of applications. Thus, there is 
now an increasing interest in also developing design procedures based on these higher-order CFD 
formulations. 6-10 If used during the design process, these higher-order CFD methods will help the 
aerodynamicist to account for the occurrence of fluid dynamic effects or phenomena which are 
not routinely predictable using potential flow methods. 

In reference 1 1 , Garabedian and McFadden described an inverse aerodynamic design procedure 
which they demonstrated using an existing FPE aerodynamics code. Their design method is based 
on a residual-correction algorithm, which we will refer to here as the GM method, and can be 
used to generate aerodynamic surfaces with prescribed surface pressure distributions. In reference 
12, Malone, et al. presented a Modified Garabedian McFadden (MGM) design algorithm that 
removed some limitations of the original GM technique. These authors applied the new MGM 
design method, also using FPE aerodynamic analysis codes as a basis, to airfoil, axisymmetric 
nacelle inlet, and 3-D nacelle inlet design problems. Later, Hazarika 13 and Sankar used a FPE 
CFD method to apply the MGM procedure to the design of blended wing-body configurations. 
In a recent effort, Malone, et al. 14 described the first use of the MGM residual-correction design 
algorithm coupled with a 2-D Navier-Stokes solution procedure. Subsequently, a similar viscous- 
flow design procedure using MGM was presented by Birckelbaw 15 , and new applications of MGM 
to multi-element airfoils using unstructured grids are under development. 16 

The objective of the present research was to develop an accurate design method for viscous, 
attached-flow, design problems which might be beyond the capability of potential-flow or Euler 
methods, even those using interactive boundary-layer theories. Because the aerodynamic designer 
normally seeks attached flow conditions, the method to be described is not expected to handle 
separated flow design problems. However, by virtue of the fact that a Navier-Stokes method forms 
the basis of the present procedure, the possible occurrence and extent of separated flow regions 
can be directly computed and noted by the designer during the design process. 

The following sections of this paper will describe the multigrid Navier-Stokes computational 
procedure, the MGM design algorithm, implementation of the design procedure, and will also 
present the results of several sample airfoil design problems. 


NAVIER-STOKES SOLUTION PROCEDURE 

The two-dimensional Navier-Stokes procedure used in the present work was originally de- 
veloped by Swanson and Turkel. 17 Their method solves the Reynolds-averaged form of the full 
Navier-Stokes equations (neglecting body forces and heat sources) on a body-fitted computational 
grid. The mathematical formulation in generalized coordinates consists of a non-dimensionalized 
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set of equations cast in a strong conservative form: 

Qt + Dc + E n = J^MooRe c (Ge H-q) 


(1) 


In Eq. (1), Q is the vector of conserved flow variables, which are themselves combinations 
of the usual primitive variables, density (p), the components of fluid velocity (u, r), and the fluid 
total energy (e). The quantity M 0 0 is the freestream Mach number and, Re c is the Reynolds 
number. The vectors D and E are the inviscid flux vectors in the f and r/ coordinate directions, 
respectively. Also, the vectors G and H are the viscous flux terms in the corresponding coordinate 
directions. The techniques used to solve Eq. (1) are given in Refs. 17, 18, and 19. Here we 
present only a brief description of the Navier-Stokes solution procedure. 

The spatial derivatives in the time-dependent Navier-Stokes equations are approximated 
with central differences. A cell-centered finite-volume technique is used to obtain the spatial 
discretization. For sufficiently smooth meshes the discretizations are second-order accurate. 
Adaptive numerical dissipation terms are appended to the resulting semidiscrete formulation. 
These terms, which are a blending of second and fourth differences, are included to provide 
shock capturing capability and to give the necessary background dissipation for convergence. In 
smooth regions of a flow field the dissipation terms are third order. The semidiscrete equations 
are integrated in time with a modified five stage explicit Runge-Kutta scheme. On the first, third, 
and fifth stages there is a weighted evaluation of the dissipation terms, which results in a good 
parabolic stability limit. The physical diffusion terms are evaluated only on the first stage and 
frozen for the remaining stages, without compromising stability. The decoupling of the temporal 
and spatial discretization makes the scheme amenable to convergence acceleration techniques, 
which are very beneficial in the computation of steady flows. 

Three techniques are employed to accelerate convergence to steady state. The first one is local 
time-stepping, where the solution at any point in the domain is advanced at the maximum time step 
allowed by stability. This results in faster signal propagation, and thus, faster convergence. The 
second technique is variable coefficient implicit residual smoothing. It can be regarded as simply 
a mathematical step applied after each Runge-Kutta stage to extend the local stability range. The 
third technique is multigrid. A multigrid method involves the application of a sequence of meshes 
to a discrete problem to accelerate convergence of the time-stepping scheme. Successively coarser 
meshes can be generated by starting with the desired fine mesh and eliminating every other mesh 
line in each coordinate direction. An equivalent fine grid problem is defined on each coarse grid. 
Appropriate operators are introduced to transfer information between the meshes. In the method 
applied here a fixed W-type cycle is used to execute the multigrid strategy. The efficiency of the 
multigrid process depends strongly upon effective high frequency damping characteristics of the 
driving scheme. Such damping behavior is provided by the five stage Runge-Kutta scheme. The 
good smoothing of the highest frequencies on the coarser meshes allows rapid removal of the low 
frequency errors in the fine grid solution. There are two additional advantages of the multigrid 
method. First, less computational effort is required on the coarser meshes. Second, information 
is propagated faster on the coarser meshes due to larger allowable time steps. 

Figure 1 presents typical computed lifts and moments for an NACA 0012 airfoil to demonstrate 
the capability of the multignd algorithm for aerodynamic analysis applications. Turbulence closure 
was obtained with the Baldwin -Lomax model. 
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MGM DESIGN PROCEDURE 

The MGM design method can be classified as a residual-correction technique, in which the 
residuals are the difference between the desired speed distribution and the computed distribution. 
Over the past decade a number of residual-correction methods have been developed, such as the 
“wavy-wall” approach of Davis . 20 The methods differ primarily in the manner in which changes 
in residual are related to changes in surface shape. The MGM algorithm itself consists of an 
auxiliary PDE that is solved for incremental changes in surface coordinates during each design 
cycle. The final aerodynamic shape is approached in a stepwise fashion through a cyclical iteration 
between the flow solver and the MGM algorithm. 

Mathematical Formulation 

The MGM auxiliary PDE is heuristic in derivation and assumes that changes in surface 
pressures are proportional to changes in airfoil surface slopes and curvatures. For two-dimensional 
flow about an airfoil configuration, the auxiliary equation is given by 

F 0 S t + F\S x t + F 2 S XIi = R (2) 

where R is the residual, defined as R = q^ — 17 /. The quantities q c and qt are the computed 
and target speed distributions, the coordinate x is the usual cartesian coordinate taken here to lie 
along the airfoil chordline, and the coefficients Fq, F\, and F 2 are constants chosen to provide 
a stable iterative process. Figure 2 shows how this auxiliary equation is typically incorporated 
into existing flow solution procedures. The computed surface velocities are normally obtained 
from partially converged numerical solutions to the flow equations under consideration at a given 
value of time, t. During the design process, as q c approaches q t , the right-hand side of Eq. (2) is 
reduced, and subsequent solutions of the auxiliary equation yield minimal changes in the airfoil 
surface coordinates. 

Next, Eq. (2) is written in terms of a correction to the airfoil coordinates, AS, by using the 
temporal derivatives and choosing At = 1, so that Eq. (2) can then be written as: 

FoAS + F,(AS), + F2(AS)„ = rt (3) 


Numerical Solution Procedure 

The auxiliary PDE is solved by writing finite-difference expressions for each term of Eq. 
(3). The computational grid used to solve this equation is the same grid used for the fluid- 
dynamic equations, which for the present Navier-Stokes solver, is an algebraically generated 
C-grid topology. Equation (3) is solved only along the airfoil surface, so that only the grid-line 
clustering in the x or streamwise direction is of importance. 

Assuming that there are a total of N computational points on the airfoil surface, Eq. (3) is 
written for each of these points, i, where 1 < i < N. A typical equation evaluated at the i th 
point on the surface is 


-t- B t A)i i -+■ CjAVi—i — R, 


(4) 
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The coefficients A„ B u and C, are evaluated by means of standard finite difference expressions, 
and AVi is the incremental change in surface coordinate, AS, at the i th computational point. 
Equation (4) is evaluated at each point, i, around the airfoil surface, leading to a system of 
equations with N unknowns, the AV, values. At each point on the aerodynamic surface, A}, is 
coupled to values at each neighboring point. The resulting algebraic equations form a tridiagonal 
system that is solved for values of AVj using the Thomas algorithm. 21 

The design cycle is completed by updating the previous surface geometry using the new 
values of A Y, as follows: 


ynew _ yold + ^ j Qr i = \ to N ( 5 ) 

Additional details of the MGM algorithm can be found in References 12 and 14. 

Trailing-Edge Crossover 

The present inverse procedure was developed to permit the design of complete airfoil surfaces, 
including the leading-edge and trailing-edge regions. However, a completely arbitrary choice for 
a target pressure distribution does not always result in a well-posed inverse design problem. 
For example, Volpe 4 has presented a technique to satisfy the three integral constraints relating 
target pressures and freestream conditions that are required to insure a well-posed problem in 
compressible flow. As a possible consequence of using unconstrained target pressures, any inverse 
procedure may produce an airfoil geometry which may exhibit trailing-edge crossover, or lead to 
other unrealistic configurations. 

Therefore an artifice is used in the present work so that the trailing edge thickness can be 
controlled and so that any tendency of the airfoil to “fish-tail” is identified. If the geometry 
is driven to a “fish-tail” configuration (trailing-edge crossover), a linear wedge is added to the 
airfoil section so that the resulting trailing-edge thickness equals a predetermined value. It has 
been demonstrated that this wedge technique can give some measure of control over the potential 
manufacturability of airfoil configurations generated by automated design procedures. 2 It should 
be noted that if the above wedging technique is required continuously during the design process, 
the original target pressures should be examined for possible modification along the lines discussed 
by Volpe 4 . A technique such as this may be used to modify these pressure distributions in order 
to rigorously provide for a well posed inverse design problem. 


RESULTS 

The MGM design procedure has been incorporated into the 2-D Navier-Stokes code described 
previously. The resulting computer program is referred to here as the MG-MGM code. In this 
section, we present three sample problems to illustrate application of the design method. Target 
pressures are obtained from a known “target geometry”, and the inverse design method is then 
used to “reproduce” the original “target” configuration. These test cases demonstrate that the 
starting geometry, or baseline configuration, used to start the design process does not have to be 
“close” in thickness or camber to the target geometry. 
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Several parameters were held constant for each of the sample problems. A “W-type” multigrid 
cycle was used throughout, together with five successive levels of grid refinement. Also, five 
“W”multigrid cycles were used between all airfoil geometry updates (ie. one design cycle). The 
computational C-grid used consisted of 321 nodes in the wrap-around, or ^-direction (33 of these 
in the wake region) and 64 nodes in the surface-normal, or ^-direction, for a total of 20,544 grid 
points. The first r/ = constant grid line was clustered to within 0.0001 chord lengths from the 
airfoil surface. Since each point on the airfoil surface is allowed to move independently, each 
can be thought of as an independent variable in the context of an optimization problem. For the 
cases presented above, there were 257 such points around the airfoil surface. 

For each case presented, a total of 160 design cycles (i.e. geometry updates) were specified. 
The program was executed on a Cray 2 and each airfoil design required approximately 16 minutes 
of CPU. Comparable Euler designs would require approximately 1 1 minutes on the same machine 
for a similarly dimensioned grid. 

Design Case No. 1 

For Case No. 1, the MG-MGM code was first used in the analysis mode to compute the 
surface pressures corresponding to an RAE 2822 airfoil at = 0.8, an angle of attack, a, equal 
to zero degrees, and Re c = 6,500,000, based on airfoil chord. 

This calculated C p distribution was then used as a target distribution for the MG-MGM code 
operated in the design mode. The baseline airfoil used to start the design was an NACA 0012 
section. As shown in Fig. 3, this airfoil is significantly different in shape from the RAE 2822 
airfoil used to produce the target pressure distribution. In this figure, as well as others depicting 
airfoil geometry, the vertical scale has been expanded. 

Figure 4 compares the design and target airfoil pressures after 40 design cycles while Fig. 
5 compares the design and target airfoil contours at this point in the design process. Figures 6 
and 7 present the corresponding comparisons for pressure and geometry after 160 design cycles. 
Figure 8 shows the results of a separate analysis computation performed after the design was 
completed. This analysis started from uniform freestream conditions (impulsive start) and used 
the grid produced by the designed airfoil contour given in Fig. 7. The comparison between design 
and target pressures is actually better than that observed during the design process. This better 
correlation exists because the pressures obtained during the design process are generated with only 
a small number of multigrid cycles on the latest computational grid. The final design corresponds 
to 160 updates to the airfoil geometry and 160 grid-generation steps. The MGM design algorithm 
itself is not computationally intensive, and because a simple algebraic grid generation scheme is 
also used in the present application, the computational overhead represents only a small fractional 
increase over that which would be required to run the original CFD method, in the analysis mode. 

Design Case No. 2 

For Case No. 2, the MG-MGM code was used in the analysis mode to compute the surface 
pressures corresponding to an NACA 0012 airfoil at M 0 0 = 0.8, an angle of attack, a = 2.0 
degrees, and Re c = 6,500,000, based on airfoil chord. 

This calculated C p distribution was again used as a target distribution for the MG-MGM code 
operated in the design mode. This time the baseline airfoil was also an NACA 0012 section. 
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However, during the design process, the freestream angle of attack was kept at a = 0.0 degrees. 
This case was used to study the possible effects that a mismatch between specified pressures and 
angle of attack might have on the design process. Figure 9 compares the baseline and target airfoil 
pressures for this case. As would be expected for a transonic flight condition, the shock wave 
locations are significantly different for the a = 2.0 targets and the a = 0.0 baseline condition. 

Figure 10 compares the design and target airfoil pressures after 40 design cycles, while Fig. 

1 1 compares the design and target airfoil contours at this point in the design process. As can be 
seen in Fig. 11, after 40 design cycles the airfoil surface has already been rotated upwards to 
adjust to the target pressure. Figure 12 presents a comparison of the geometry after 160 design 
cycles. As in the previous case, a separate analysis run was performed to verify the airfoil design. 
Figure 13 shows the results of the separate analysis computation performed after the design was 
completed. This analysis started from uniform freestream conditions (impulsive start) and used 
the grid produced by the designed airfoil contour given in Fig. 12. Finally Fig. 14 shows a 
plot of the average A q 2 versus multigrid work for the 800 multigrid cycles. This quantity drops 
approximately two orders of magnitude during the design process and is used to monitor the 
progress of the design algorithm. 

Design Case No. 3 

The final example problem, design Case No. 3, was chosen to demonstrate that large geometric 
changes can be achieved with the MGM design algorithm. For this application, the target pressures 
corresponded to an NACA 0012 airfoil at Moo= 0.8, angle of attack, a = 0.0 degrees, and Re c 
= 6,500,000, based on airfoil chord. The baseline configuration used was an NACA 0006 airfoil. 
A comparison of the target and final design airfoil shapes is shown in Fig. 15. A comparison 
of the target pressures, and those obtained from a separate analysis (impulsive start) of the final 
design configuration are shown in Fig. 16. In this example, an airfoil design was successfully 
accomplished which required a 100% increase in airfoil thickness over that of the baseline airfoil 
shape. 


CONCLUDING REMARKS 

The MGM design procedure has been incorporated into an existing multigrid Navier-Stokes 
code. The computational efficiency of the method indicates that it is a viable tool for the 
design process. The actual computational effort of this design method depends, of course, on 
the complexity of the target pressure distributions chosen. Normally, aerodynamicists would seek 
to eliminate shockwaves due to the impact of wave drag on performance. Previous experience 
with the MGM algorithm 14 indicates that shock-free design applications require about 50% less 
computational effort than for flows with shockwaves present. The transonic flow cases shown 
here were picked, in part, to demonstrate the design algorithm’s robustness and ability to respond 
correctly to shockwaves in the flowfield. This feature is important because regions of sonic flow 
may be created locally near regions of high airfoil curvature even at relatively low freestream 
Mach numbers. 

Because of the computer resource requirements, any Navier-Stokes based design method would 
likely be used in combination with other, lower-cost design methods. For example, an initial 
airfoil shape designed with a FPE method may prove to be an excellent starting configuration for 
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a higher-order design approach. Used in this manner, the present Navier-Stokes inverse design 
method should then be able to account for viscous flowfield phenomena that may not be detected 
or predicted accurately enough by other methods based on FPE or Euler solution procedures. 
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(a) Lift curve (b) Drag polar 

Figure 1. Variation of lift and drag coefficients for NACA 0012 airfoil 
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Figure 2. Implementation of the MGM Inverse Design Algorithm 



Figure 3. Comparison of target and baseline 
airfoil contours for case no. 1 



Figure 4. Comparison of target and inverse 
airfoil pressures for case no. 1 
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Figure 5. Comparison of target and inverse Figure 6. Comparison of target and inverse 

airfoil contours for case no. 1 airfoil pressures for case no. 1 
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igure 1 1 . Comparison of target and inverse 
airfoil contours for case no. 2 


Figure 12. Comparison of target and inverse 
airfoil contours for case no. 2 
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Figure 13. Comparison of target and final 
airfoil pressures for case no. 2 
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Figure 14. Average Aq2 versus multigrid 
work performed for case no. 2 
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Figure 15. Comparison of target and inverse 
airfoil contours for case no. 3 


Figure 16. Comparison of target and final 
airfoil pressures for case no. 3 








